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Abstract.  The  usual  way  of  computing  partial  correlations  is  based  on  the  formation  of  the 
covariance  matrix,  that  amounts  to  squaring  the  data  matrix,  thus  inviting  a  potential  loss  of 
numerical  accuracy.  This  paper  recommends  the  determination  of  partial  correlations  from  the 
data  matrix:  the  QR  decomposition  of  the  data  matrix  is  computed  and  plane  rotations  are  applied 
to  the  resulting  upper  triangular  matrix,  which  is  the  Cholesky  factor  of  the  covariance  matrix. 
We-show  that  if  rotations  are  applied  to  the  triangular  matrix  so  as  to  leave  the  number  of  its 
zero  entries  invariant,  the  sines  of  the  rotation  angles  are  partial  correlations.  Different  ways  of 
organizing  the  computations  are  presented  for  extracting  any  set  of  partial  correlations. 
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1.  Introduction 


Classical  texts  [l]  on  multivariate  statistics  suggest  the  computation  of  partial  correlations  by 
first  forming  the  empirical  covariance  matrix  j^tj(A  —  A)T{A  -  A),  where  A  is  the  m  x  n  data 
matrix,  whose  ith  column  is  associated  with  random  variable  A4,  and  A  =  ^eeT A  is  the  empirical 
mean  matrix  ( e  is  the  m  x  1  vector  of  ones).  Substantial  loss  of  numerical  accuracy  is  incurred 
by  squaring  the  data  matrix  thus  resulting  in  errors  in  the  computed  partial  correlations  (one  can 
easily  construct  examples  where  a  full-rank  matrix  A— A  leads  to  a  numerically  indefinite  covariance 
matrix).  This  loss  of  accuracy  is  inherent  in  the  use  of  the  covariance  matrix  and  independent  of 
the  formulas  and  recursions  employed  to  compute  the  partial  correlations. 

These  shortcomings  may  be  avoided  with  a  method  working  directly  on  the  data  matrix  and, 
in  addition,  employing  orthogonal  transformations.  Our  approach  consists  of  two  steps:  after  the 
QR  decomposition  of  the  matrix  A  —  A  the  resulting  upper  triangular  matrix  U  is  transformed  to 
lower  triangular  form  L  via  plane  rotations.  The  rotations  are  executed  in  a  specific  order  which 
exploits  the  zero  structure  of  the  upper  triangular  matrix,  and  the  values  of  their  sines  constitute 
the  partial  correlations  between  variables  A,  and  Aj,  for  »  <  j,  holding  variables  A,+i  .  ..Ay_i 
fixed. 

By  effecting  the  rotations  on  only  a  submatrix  of  U,  the  partial  correlations  between  A,  and 
Aj  where  variables  At . . .  A,_i  and  A,>i  . . .  Ay_j  are  held  fixed  are  efficiently  computed,  without 
having  to  reorder  the  columns  of  either  the  data  matrix  A  or  the  upper  triangular  matrix  U.  In 
general,  we  will  present  ways  of  organizing  the  computations  so  as  to  determine  any  set  of  partial 
correlations  while  keeping  arbitrary  collections  of  variables  fixed . 


2.  Partial  Correlation  Coefficients 


The  column  vector  of  m  observed  values  a*,,  1  <  k  <  m,  of  a  real  random  variable  A,  is 
denoted  by 


o,  = 


aw 

an 


The  centered,  or  zero-mean,  data  vector  is  a*  =  ay  —  a*,  where  the  barred  quantity  denotes  the 
mean  vector 


a,  =  — e(e7'oi), 
m 


and  e  is  the  column  vector  of  m  ones. 


The  empirical  correlation  coefficient  pi,  of  two  random  variables  Ai  and  Ay  is  defined  as  the 
cosine  of  the  angle  9H  between  the  centered  data  vectors  a,  and  ay : 

Pi)  =  COS  0,y  = 

The  correlation  between  two  variables  A\  and  Aj  can  arise,  in  part,  from  the  fact  that  both 
Ai  and  Aj  show  a  correlation  with  a  third  variable  A*.  The  ‘partial  correlation  between  A,  and  Aj 
given  Ak *  then  represents  the  correlation  between  A,  and  Ay  after  the  dependence  on  At  has  been 
removed.  Formally,  the  empirical  partial  correlation  coefficient  p^-  between  variables  A,  and  Ay 
given  ( conditioned  with  respect  to)  variable  A*  is  defined  to  be  the  cosine  of  the  angle  between 
a*  and  a*  where 

af  =  ot i-  (afat)(aj’at)_1at,  a£  =  a,  -  (aj'at)(aj'at)"1ajt) 

are  the  respective  projections  of  a,  and  ay  onto  the  subspace  orthogonal  to  the  vector  otk  (note 
that  superscripts  here  denote  conditioning  rather  than  powers).  Substituting  this  into  the  formula 
for  the  cosine 

.  («? >r«; 

yields 

pk  _  aiai  ~ 

yjajai  -  (atTa* )  (a£a* )  “ 1  (a[ a,) ^ ay  -  (aj'ai)(aja*)-1(a*  ay) 

Note  that  all  quantities  in  the  expression  for  p,*  are  of  the  form  atray.  This  means  that  in  the 
general  case  of  n  variables  Ai  with  m  observed  data  values  each,  the  partial  correlations  can  be 
computed  from  the  elements  of  the  n  x  n  empirical  covariance  matrix 

B  =  ( bij )  =  (A  -  A)T (A  -  A), 

where  A  -  A  is  the  centered  m  x  n  data  matrix,  and  the  data  matrix  A  and  its  mean  matrix  A  are 
defined  by 

A=[oi  ...  an  j ,  A  =  —  e{eT  A) 

m 

(the  empirical  covariance  matrix  is  usually  defined  as  (m  -  l)-1  B;  since  partial  correlations  are 
normalized  quantities,  independent  of  the  scaling  by  (m— 1)_1 ,  we  shall  use  here  the  more  convenient 
unsealed  expression).  Denoting  by 


a;  a. 


bij  =  bij  -  bikb^bkj  =  (af ay  -  (aj'at)(aj'at)  ‘(ajay))  =  (a‘)7’a* 


the  elements  of  the  Schur  complement  in  B  with  respect  to  [2,  3]  one  has 


nn 

In  general,  the  conditioning  may  occur  with  respect  to  more  than  cne  variable,  for  instance, 
with  respect  to  Ai,  A*  and  Am  or  with  respect  to  a  sequence  A*  . . .  A*+j.  In  that  case  the  involved 
vectors  a,-  and  ay  are  projected  onto  a  subspace  orthogonal  to  the  subspace  spanned  by  a*,  at 
and  am  or  by  a*  . . .  a*+j ,  respectively.  Denoting  by  bk'l,m  and  6kjk+1  the  elements  of  the  respective 
Schur  complements  of 


■ 

bu 

bkm 

bkk  ■ 

bk,k+l 

bik 

bu 

blm 

and 

-  bmk 

bml 

bmm  • 

■  bk+iik+i . 

in  B,  the  partial  correlation  between  Ai  and  Aj  given  At,  Ai  and  Am  is 

_  °ij 

and  the  partial  correlation  between  A,  and  Ay  given  A*,  A*+1, . . . ,  Ai+(  is 

bk:.k+l 

s  _ 


Our  notation  automatically  incorporates  the  so-called  quotient  property  for  Schur  complements 
[2],  which  essentially  states  that  the  effect  of  conditioning  with  respect  to  variables  belonging  to  a 
set  S  can  be  accomplished  by  first  conditioning  with  respect  to  variables  that  belong  to  a  subset  of 
Sj  of  5  followed  by  conditioning  with  respect  to  the  remaining  variables  in  S  -  Si ,  the  complement 
of  Si  in  S.  The  quotient  property  for  Schur  complements  yields  readily  recursive  formulas  for  the 
computation  of  pky,,m  or  pk^k+1  and  these  formulas  are  the  ones  generally  used  to  compute  the  partial 
correlation  coefficients  [1,  3j.  An  inconvenience  with  such  formulas,  that  rely  on  the  computation 
of  Schur  complements  in  the  covariance  matrix,  is  that  construction  of  the  covariance  matrix  itself 
implies  squaring  up  the  data,  A  -  A,  and  thus  a  doubling  of  the  dynamic  range  and  potential  loss 
of  accuracy  as  the  subsequent  example  shows. 

Example.  Suppose  that  the  following  zero-mean  data  matrix  is  given, 

■  1  -1  c 

.  i  -1  1  -c 

A-A=-= 

V2  e  e  1 

,-e  —  c  -1. 
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where  e  is  non*zero,  and  that  the  partial  correlation 


Pis  = 


is  to  be  determined.  The  corresponding  covariance  matrix  is 


-  1+e2 

-1+e2 

2e  - 

B  = 

-1  +  e2 

1  +  e2 

0 

.  2c 

0 

1  +  e2. 

and  in  exact  arithmetic,  one  has 


so  that 


l2  _  n  l2  _  ,  .  r2  +  .2  .  ,  .  2 

"13  ~  *e,  ®11  —  1  +  e  j  ~  f2 - »  O33  —  1  +  f 


2  f  1  C  >  0 

Pis  =  «gn(<)  =  S 

l-l  e  <  0 


However,  in  finite  precision  floating  point  arithmetic  and  with  c  chosen  to  be  sufficiently  small  (e.g. 
e  is  the  largest  number  so  that  the  computed  fl(l  ±  4c3)  =  1),  the  computed  quantities  turn  out  to 
be 

•  1  -1  2C 

fl(B)=  -1  1  0 

.2c  0  1  . 

and 

^13  =  2e,  b\  j  =  0,  633  =  1 

so  that  p\ 3  is  not  a  finite  number. 

The  next  section  introduces  a  numerical  method  that  achieves  much  higher  accuracy  by  working 
directly  on  the  data  matrix  A  —  A. 

In  order  to  avoid  squaring  the  data  matrix  one  may  try  to  work  with  the  Cholesky  factor  of  B. 
Indeed  the  n  x  n  upper  triangular  Cholesky  factor  U  can  be  obtained  without  squaring  from  the 
QR  factorization  of  the  scaled  centered  data  matrix  since,  if  A  -  A  is  decomposed  into  the  product 
QU  where  Q  has  orthogonal  columns, 

B  =  (A  -  A)t(A  -  A)  =  (A  -  A)tQQt(A  -  A)  =  UTU. 

It  is  known  that  the  elements  of  the  Cholesky  factor  can  be  represented  in  terms  of  elements 
of  Schur  complements  with  respect  to  the  leading  principal  submatrices  of  B : 

4 


Lemma  2.1.  The  aoa-zero  elements  of  the  Cholesky  factor  U  of  B  are  of  the  form 

j  >  i. 

Proof.  Let  U  be  an  upper  triangular  matrix  with  elements  u,y  =  ,  j  >  ».  The 

(i,j)th  element,  j  >  t,  of  the  symmetric  matrix  UTU  is 

k=l  k=l  k=l 

= -  (»<,  -  i,i(tn)-‘t.;) + E‘«  "‘(‘if'1)''6!/ ■■ = *. -  +e*« 

i=2  i=2 

= - »?,«,)-%) + ZW-'vit-'rxr' = *«  -  V + ilC-'ttl-'r'iir 

k=i  k=3 

-•••=<«-  (‘ir1  -  w-'pj-'rvt')  - 

where  the  telescoping  of  the  sum  is  achieved  by  making  use  of  the  quotient  property  of  Schur 
complements.  Hence  CPv  =  H,  and  the  uniqueness  of  the  Cholesky  factor  implies  U  =  U. 


The  partial  correlation  between  A^  and  Ay  given  the  intermediate  variables  Aj,  . . . ,  A,_i  can 
be  expressed  as 

»!r'  =  = 

Thus  it  seems  that  the  partial  correlations  may  be  computed  as  simple  functions  of  the  elements 
of  the  Cholesky  factor.  Yet  unfortunately  the  desired  quantity  (bj? -1)-1/2  differs  from  u^1  = 
Moreover,  it  is  hard  to  see  how  to  determine  the  quantities  (fcV* without  the 
use  of  squaring  operations. 

The  next  section  introduces  a  numerical  method  that  gets  around  this  difficulty  by  applying 
plane  rotations  to  the  columns  of  the  Cholesky  factor. 

3.  New  Algorithm 

Prom  the  previous  section  it  is  clear  that  one  must  think  of  more  subtle  means  to  employ  the 
Cholesky  factor:  our  method  determines  partial  correlations  as  cosines  evaluated  through  inner 
products  but  instead  as  sines  of  rotations  that  zero  out  components  of  certain  column  vectors.  The 
key  idea  for  the  new  algorithm  is  based  on  the  fact  that  the  data  vectors  are  initially  represented 
by  a  triangular  matrix,  the  Cholesky  factor;  and  that  the  partial  correlations  may  be  computed  by 


kl'Ll'Ll'LlW  V  :i 


ir/2  -  0i2 


Figure  1:  Angles  in  the  2x2  Example. 

applying  plane  rotations  in  a  particular  order  to  the  columns  of  the  Cholesky  factor.  To  see  that 
consider  a  simple  2x2  example. 

Let 

Uii  Ui2 

u  = 

0  ti22 

be  the  upper  triangular  factor  (with  positive  diagonal  elements)  in  the  QR  decomposition  of  a 
m  x  2  matrix  A  -  A.  The  (partial)  correlation  pu  between  Aj  and  A2  is  the  cosine  of  the  angle  0l2 
between  the  two  columns  of  U .  Because  of  the  triangular  structure  of  U  its  first  column,  [  un  0  ]T, 
is  a  positive  multiple  of  the  first  canonical  vector  ei  =  [  1  0  ]T  while  its  second  column  is  a  linear 

combination  of  ex  and  the  second  canonical  vector  e2  =  [0  1  ]T. 

The  columns  of  the  matrix  U  may  be  rotated  in  such  a  way  that  the  second  column  becomes 
a  positive  multiple  of  e2  thereby  turning  the  first  into  a  linear  combination  of  e\  and  e2: 


L  =  QU  = 


*  n  0 
*21  *22 


Suppose  the  angle  between  t\  and  e2,  denoted  by  ^(ei,e2),  is  +jr/ 2  then  the  angle  between  the  two 
columns  can  be  defined  as  0i2  =  The  fact  that  the  first  column  is  a  positive  multiple  of 

ei  implies  ^(ei.uj)  =  0i2-  To  turn  the  second  column  into  a  positive  multiple  of  e2  requires  that 
all  columns  of  U  be  rotated  by  the  angle 

Z(u2,e2)  =  ^(ei,e2)  -  /.{ei,u2)  =  jt/2  -  9U, 

see  Figure  1.  Hence  the  angle  between  the  two  columns  of  U  is  preserved  under  the  rotation,  and 


$ 


i 


the  angle  of  such  a  rotation 

c  — 8 

e  = 

8  C 

completes  £«  to  a  right  angle:  c  =  cos  (jt/2  -  #12)  and 

8  =  sin  (jt/2  -  #12)  =  cos On  =  Pn- 

Consequently,  the  desired  (partial)  correlation  is  the  sine  of  the  rotation  0. 

The  above  suggests  that,  in  general,  certain  partial  correlations  may  be  computed  from  the 
plane  rotations  that  transform  the  upper  Cholesky  factor  to  the  lower  Cholesky  factor.  The  tri¬ 
angular  zero-structure  of  U  makes  it  possible  to  rotate  columns  in  a  manner  illustrated  above  and 
determine  a  partial  correlation  from  the  sine  of  a  rotation, 
brief  look  at  the  3x3  case 

'Uji  U12  uis' 

U  =  0  1*22  U23 

.  0  0  U33 . 

illustrates  the  above. 

At  first,  because  the  second  column  has  only  one  more  non-zero  element  than  the  first,  the 
columns  of  U  can  be  rotated  in  the  (ei,e2)-plane  so  as  to  make  the  second  column  co- linear  with 
«2, 

r  C12  -*12  01  run  uj2  ujs  1  r*  0  *  1 


‘  * 

0 

*  “ 

= 

* 

* 

♦ 

.0 

0 

U33- 

®12  C12  0  0  U22  U2S  =  *  *  *  , 

.0  0  1.  .  0  0  U33 .  .00  U33 . 

and  P12  =  S12.  Here,  ‘co-linear’  is  used  to  mean  ‘a  positive  multiple  oP  and  *  denotes  terms  that 
are  non-zero  in  general. 

Next,  to  achieve  conditioning  of  Ai  and  A3  with  respect  to  A2,  the  first  and  third  columns 
are  projected  onto  the  subspace  orthogonal  to  the  second  column.  Due  to  the  triangular  structure 
of  U  and  the  effect  of  the  previous  rotation  the  second  column  is  co-linear  to  C2>  and  the  subspace 
orthogonal  to  it  is  just  the  plane  (ei,  es).  The  partial  correlation  p\3  can  then  be  determined  from 
that  rotation  that  makes  co-linear  with  63  the  projection  of  the  third  column  onto  (ei,es).  Since 
this  rotation  takes  place  in  a  subspace  orthogonal  to  the  second  column  it  does  not  affect  the  second 
column,  and  the  zero  element  introduced  by  the  previous  rotation  is  preserved: 


0  -«?, 


0  1 
-®13  0 


■*  0  *  • 

*  *  *  = 

.0  0  U33 . 
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and  />ig  =  <fg.  Note  that  another  non-zero  element  is  introduced  in  the  first  column. 

Again,  because  of  the  triangular  structure  of  U  and  the  effect  of  the  second  rotation  the  zero- 
structure  of  the  second  and  third  columns  is  the  same  save  for  one  element,  the  second  column 
is  co-linear  with  e j  while  the  third  is  a  linear  combination  of  e 2  and  e$.  Thus  the  columns  of 
the  matrix  can  be  rotated  to  yield  P23  by  applying  a  rotation  that  makes  the  whole  third  column 
co-linear  with  es,  and  turns  the  second  column  into  a  linear  combination  of  e j  and  ej. 


■10  0  - 

1 - 

O 

O 

w< 

*■4 

'In  0  0  - 

0  C2S  -S2S 

*  *  * 

= 

hi  hi  0 

.0  «23  C 23  - 

- 1 

* 

O 

* 

_ 1 

•  hi  hi  hs  • 

and  p 23  =  «23- 

Theorem  3.1.  If  the  elements  in  the  Cholesky  factor  U  of  the  covariance  matrix  B  are  eliminated 
in  the  order 

■*  1  2  ...  n  -  1  ' 

*  n  2n  -  3 

*  *  > 

*  n(n  -  l)/2 
* 

that  is,  proceeding  row  after  row  from  top  to  bottom,  and  within  each  row  from  left  to  right,  then  the 
sine  of  the  rotation  that  eliminates  element  ( i,j ),  j  >  i,  is  equal  to  the  partial  correlation 

Proof.  The  proof  proceeds  by  induction. 

The  induction  basis  comprises  the  computation  of  partial  correlations  between  A\  and  all  other 
variables.  To  start  with,  the  matrix  U  is  of  the  form 

'*  *  *  *  ...  *” 

*  *  *  ...  * 

*  *  ...  * 

*  ...  * 


From  the  2x2  case  one  can  see  that  elimination  of  element  (1,2)  in  the  upper  triangular  matrix  U 
by  a  rotation  in  plane  (ej,C2)  provides  pw  The  second  column  of  the  resulting  matrix  becomes 
co-linear  to  e2  while  the  first  column  becomes  a  linear  combination  of  ej  and  e2-  Hence,  there  is  a 


new  non-zero  element  in  the  first  column  and  a  zero  has  been  introduced  in  the  first  row: 


“  *  *  *  ...  * " 

*  *  *  *  ...  * 

*  *  ...  * 

*  ...  * 

* . 

The  3x3  case  showed  that  the  correlation  p\z  between  A\  and  A$  given  Ai  could  be  computed  by 
rotating  the  first  and  third  column  and  thereby  introducing  a  non-zero  element  in  position  (3, 1) 
and  a  zero  in  position  (1,3): 


*  *  *  *  ...  * 

*  **...* 

*  ...  * 

* . 

Continuing  this  argument,  the  partial  correlation  p1;-  between  A\  and  Aj,  given  Ai>  . . Aj~\  is 
computed  by  peforming  a  rotation  in  plane  (ei,e;)  thereby  creating  a  zero  element  in  position  (1,  j). 
Thus,  once  all  correlations  involving  A\  have  been  computed  the  first  column  of  the  matrix  has 
totally  filled  in,  and  the  first  row  is  zero  except  for  the  first  element: 


*  *  *  * 

*  *  * 

♦  * 

* 


♦ 

* 

* 


Assume  that  the  partial  correlations  pkkfll~l 


have  already  been  computed  for  k  <  i  and 


Vii 

vi+ 1,1+1 

V1+1,J-1  vi+l  J 

Vi-U- 1  Vi-U 

viJ 

where  Lq  is  lower  triangular  and  Uq  is  upper  triangular. 

By  induction  hypothesis  the  entire  lower  triangular  part  of  the  leading  i  -  1  columns  is  non¬ 
zero,  and  the  ith  column  has  j  —  2  non-zeros  in  its  lower  triangular  part  due  to  the  computation  of 
Pi,i+i,  ■  •  ■  >  P*ij Ii~2-  In  order  to  compute  the  next  correlation  the  corresponding  columns 

Vi,  . ..,  Vj  of  the  current  matrix  must  be  projected  onto  a  subspace  orthogonal  to  the  subspace 
spanned  by  Aj+t,  . . . ,  Ay_i.  Due  to  the  initial  ‘nesting’  of  the  column  subspaces  (i.e.  the  original 
upper  triangular  structure  of  U)  the  trailing  components  j,  . . . ,  n  of  vt+1,  . . . ,  are  zero;  and 
due  to  the  rotations  performed  in  order  to  retrieve  previous  partial  correlations  (i.e.  the  appearing 
lower  triangular  structure  of  L)  the  leading  components  1,  ...,  i  of  v,-+1,  ...,  vy_i  are  zero. 
Hence  the  subspace  spanned  by  A,+1,  . . . ,  Ay_i  is  the  space  spanned  by  e,+i,  . . . ,  ey_i,  and  the 
space  orthogonal  to  it  is  the  space  spanned  by  ej,  ...,  e*,  e;,  ...,  e„.  Similarly,  components 
1,  . . . ,  »  —  l,  j,  . . . ,  n  of  and  components  1,  . . . ,  t  —  1,  j  +  1,  . . . ,  n  of  Vj  are  zero;  and  the 
projections  of  Vj  and  onto  ex,  . .  ,  e,,  ey,  . . . ,  e„  are  respectively  co-linear  to  e,  and  a  linear 
combination  of  e,  and  ey.  Thus,  p\*x  J_1  is  obtained  by  applying  the  rotation  in  plane  (e,-,  ey)  that 
makes  the  projection  of  vy  co-linear  with  e*;  p’*1 *_1  is  the  sine  of  that  rotation.  After  the  rotation 
the  matrix  has  the  form 
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Remark.  If  the  matrix  L  is  not  needed  about  half  of  the  arithmetic  operations  can  be  saved  by 
applying  the  rotations  merely  to  the  trailing  principal  submatrix  of  interest. 

Returning  to  the  example  of  the  previous  section  it  becomes  clear  that  the  new  method  can 
avoid  the  loss  of  accuracy  associated  with  the  explicit  formation  of  the  covariance  matrix. 

Example.  Performing  a  QR  decomposition  of  the  matrix  A  -  A  yields  4x3  matrix  Q  with 
orthonormal  columns  and  a  3  x  3  upper  triangular  factor 

'  I  -f  c2  - 1  +  c2  2c 

U  =  0  2'c'  “*“(e)(l  - c2) 

Loo  0 


_1 _ 
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in  exact  arithmetic  and 

■1  -I  2c  ' 
fl(t/)  =  0  2\c\  sign(t) 

.0  0  0  . 

in  finite  precison  arithmetic  with  the  same  choice  of  c  as  before.  The  first  rotation  with  fl(si2)  =  -1 
and  fl(ci2)  =  2\c\  yields 

r2|f|  0  sign(e)‘ 
fl(tf')=  -1  1  0 

.0  0  0  . 

and  the  second  rotation  gives  cj3  =  0  and  sj3  =  sign(e)  so  that  fl(pj3)  =  sign(e)  indicates  a  linear 
dependence  between  the  three  columns  of  A  -  A  as  in  the  true  computation. 

4.  Computation  of  Arbitrary  Partial  Correlations 

Subject  to  a  certain  inital  ordering  of  the  random  variables  A\ ,  . . . ,  A„  our  algorithm  computes 
the  partial  correlations  between  A,  and  A,  given  A*+i,  ...,  Ay_i  by  completely  reducing 

the  upper  triangular  matrix  U  to  a  lower  triangular  matrix  L. 

Other  partial  correlations  may  be  computed  by  performing  only  a  partial  reduction.  For 
instance,  consider  the  following  6x6  example 
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The  leading  three  columns  of  U  span  the  subspace  of  A\ ,  Aj  and  As,  and  this  is  equal  to  the  space 
spanned  by  the  first  three  canonical  vectors  t\,  e2  and  ej  due  to  the  triangular  structure  of  U .  The 
space  orthogonal  to  it  is  the  one  spanned  by  £4,  e&,  eg  and  is,  because  of  the  triangular  structure, 
equal  to  that  of  columns  4  through  6  of  U  with  components  1  to  3  set  to  zero.  This  means  that 
the  correlation  between  A4  and  As,  given  A\,  As  and  As,  can  be  computed  by  a  rotation  of  U 
in  plane  (e4,es).  The  resulting  matrix  has  a  new  zero  in  column  five  and  a  fill-in  in  column  four: 
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The  next  correlation  that  can  be  computed  is  pig'5  with  a  rotation  in  plane  (e4,e6),  the  subspace 
orthogonal  to  Aj,  A2,  A3  and  As: 
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The  last  correlation  p\$  is  determined  by  completing  the  transformation  of  the  3x3  trailing 
principal  submatrix  to  lower  triangular  form: 
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In  general,  the  correlation  p^a',+l^  1  for  1  >  a  and  i  <  j  <  n  can  be  determined  be  preserving 
the  leading  a  rows  and  columns  of  U  and  transforming  the  trailing  principal  submatrix  of  order  n-cx 
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If  it  is  known  in  advance  which  partial  correlations  are  to  be  determined  then  the  columns 
of  the  m  x  n  data  matrix  may  be  ordered  so  as  to  minimize  the  number  of  arithmetic  operations 
succeeding  the  computation  of  the  Cholesky  factor. 

For  instance,  a  lower  bound  on  the  number  of  arithmetic  operations  in  the  computation  of 
pfj,  where  5  is  a  subset  of  k  >  0  numbers  in  1  ...  n  not  containing  i  and  j,  is  0(n  —  k)  since 
our  method  requires  at  least  one  rotation  to  compute  a  partial  correlation  and  the  dimension  of 
the  space  involved  is  n  -  k.  This  lower  bound  is  attained  by  ordering  the  columns  so  that  the 
set  5  represents  the  leading  k  columns  of  the  data  matrix  followed  by  columns  A<  and  Aj.  The 
correlation  p *•  can  then  be  determined  by  one  rotation  in  the  plane  (efc+i,e*+2)  that,  due  to  the 
triangular  structure  of  the  Cholesky  factor,  involves  0(n  -  k)  non-zero  element  pairs. 

Not  only  the  ordering  of  the  columns  is  important  but  also  the  sequence  in  which  particular 
correlations  are  computed.  Consider  the  computation  of  a  partial  correlation  between  two  variables 
Ai  and  Aj  with  successively  more  variables  fixed:  pfl,  . . .,  pff ,  where  Si  C  ...  C  Sk  and  t,  j  g  Sk¬ 
it  seems  that  the  following  order  of  rotations  constitutes  the  simplest  way  of  determining  the  above 
correlations.  It  is  illustrated  by  means  of  a  5  X  5  example  for  the  computation  of  ^>12,  p*2,  Pi  *,  and 
Pi At  first  the  columns  of  the  data  matrix  are  ordered  so  that  i  and  j  represent  the  first  two 
columns  followed  by  the  columns  of  Si,  the  columns  of  Sj  —  Si,  the  columns  of  S3  —  S3  —  Si,  etc.  In 
the  example  this  amounts  to  the  ‘natural’  ordering  A\  . . .  As  of  the  variables.  The  first  correlation 
P12  can  now  be  computed  with  one  rotation  from  the  Cholesky  factor  U.  To  compute  p\2  columns 
two  and  three  of  U  are  exchanged  and  a  rotation  in  plane  (e2,es)  results  in  the  Cholesky  factor  U' 
corresponding  to  the  data  matrix  with  variables  in  the  order  Ai,  A3,  A2,  A4,  A5.  Since  A3  is 
situated  between  Ai  and  A2  two  rotations  suffice  for  the  computation  of  p\2.  The  effect  of  these 
steps  on  the  matrix  is  depicted  below: 
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Similarly,  exchanging  columns  three  and  four  of  U',  performing  a  rotation  in  plane  (es,C4)  to  get 
the  Cholesky  factor  U"  of  the  data  matrix  corresponding  to  the  ordering  Aj,  As,  A4,  A2,  A5,  and 
performing  three  more  rotations  on  U"  results  in  the  extraction  of  In  general,  if  the  sets  Si 
differ  by  more  than  one  index,  more  columns  of  the  Cholesky  factor  must  be  exchanged  to  ensure 
that  all  fixed  variables  are  situated  between  A,  and  A,. 

As  for  arbitrary  sequences  of  partial  correlations,  the  determination  of  the  column  ordering  of 
the  data  matrix  as  well  as  the  computation  sequence  of  the  partial  correlations  so  as  to  minimize 
the  number  of  arithmetic  operations  seems  to  be  an  NP-complete  problem.  The  use  of  heuristics, 
such  as  the  following  greedy  approach,  might  lead  to  acceptable  operation  counts:  the  random 
variables  are  ordered  so  that  as  many  partial  correlations  as  possible  can  be  determined  from  the 
resulting  Cholesky  factor.  Repeatedly,  the  columns  of  the  Cholesky  factor  are  then  re-ordered 
according  to  the  same  strategy,  the  matrix  returned  to  upper  triangular  form,  and  appropriate 
rotations  performed  until  all  correlations  have  been  computed. 

5.  An  Open  Problem 

Adding  a  row  aT  to  the  m  x  n  data  matrix  A  results  in  a  rank-two  update  to  the  Cholesky 
factor  U  of  A  -  A.  Suppose  the  QR-factorization  of  A  -  A,  A  -  A  =  QU  where  Q  is  m  x  n  with 
orthogonal  columns  is  available.  At  first  the  rank-one  n  x  n  matrix  QTe( ^ 

which  can  be  computed  in  0(mn)  operations,  is  added  to  U  and  then  the  row  aT  -  ^j(er  A-f-  aT) 
is  appended  as  the  (n  +  l)st  row.  The  so  augmented  Cholesky  factor  can  be  reverted  to  upper 
triangular  form  U'  by  means  of  rotations  in  0(n2)  operations.  The  quantities  of  interest,  the  partial 
correlations  of  the  updated  matrix,  can  be  computed  from  U' .  However,  instead  of  starting  all  over 
from  U',  a  more  effective  approach  could  be  to  use  the  augmented  Cholesky  factor  as  starting  point 
in  order  to  ‘update’  the  partial  correlations  computed  from  U .  A  similar  problem  arises  after  the 
deletion  of  a  row  from  the  data  matrix. 
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